#ifndef MATRIX_H
#define MATRIX_H

#include <vector>

#include "Vector/vector.03.h"
#include "GadgetParticleType/gpt.01.h"

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>



class Matrix { 
	
	public: 
	
	static const double pi = 3.1415926535897932384;		
	
	dVector X, Y, Z;	
	
	// CALCULATION OF INERTIA TENSOR
 	void calc(uint N, P *p) { 	
		
		X=Y=Z=dVector(0,0,0);
		for (int i=0; i<N; i++)
			add_one(p[i].Mass, *(fVector*)p[i].Pos);
	}
	
	void calc(uint N, float *m, Vector<float> *r) { 	
		
		X=Y=Z=dVector(0,0,0);
		for (int i=0; i<N; i++) add_one(m[i], r[i]);
		
	}
	
	template <typename T1,typename T2> 
	void calc(vector<T1> &m, vector<T2> &r) { 	
		
		X=Y=Z=dVector(0,0,0);
		for (int i=0; i<r.size(); i++) add_one(m[i], r[i]);
		
	}
	
	template <typename T> 
	void calc(vector<T> &p) { 	
		
		X=Y=Z=dVector(0,0,0);
		for (int i=0; i<p.size(); i++) add_one(p[i].m, p[i].R);
		
	}	
	
	void calc(uint N, float m, Vector<float> *r) { 	
		
		X=Y=Z=dVector(0,0,0);
		for (int i=0; i<N; i++) add_one(m, r[i]);
		
	}

	template <typename T1,typename T2> 
	void add_one(T1 &m, Vector<T2> &r) { 	
		
		X +=m*Vector<T2>( SQR(r.y) + SQR(r.z),  -r.x*r.y, r.x*r.z );
		Y +=m*Vector<T2>( -r.y*r.x, SQR(r.x) + SQR(r.z),  -r.y*r.z );
		Z +=m*Vector<T2>( -r.z*r.x, -r.z*r.y,  SQR(r.x) + SQR(r.y) );
		
	}
	
	
	
	
	
	dVector A, B, C;	
	
	void calcEigenVectors(){
		
		uint N = 1500000;
		double Rmaxmod=0, Rminmod=1e99;
		dVector Rmax, Rmin;
		for (uint i=0; i<N; i++) {
				
			double x = 1-2*(double(i)/N);
			double cosa1 = x;
			
			double a1 = acos(x);
				
			double a2 = pi*((double)rand()/RAND_MAX);	
			double sina1 = sin(a1);
				
			dVector R = cosa1*X + sina1*(cos(a2)*Y + sin(a2)*Z);
			double Rmod =R.module();
				
			if (Rmod > Rmaxmod) { Rmax = R; Rmaxmod = Rmod;}
			if (Rmod < Rminmod) { Rmin = R; Rminmod = Rmod; }
				
		}
		
		A = (Rmax).ort();
		B = (Rmax%Rmin).ort();
		C = (Rmin).ort();
	}
	
	void calcEigenVectorsGSL(){
		
		double data[] = {	X.x, X.y, X.z, 
								Y.x, Y.y, Y.z,
								Z.x, Z.y, Z.z };
								
		gsl_matrix_view m = gsl_matrix_view_array (data, 3, 3);
		
		gsl_vector *eval = gsl_vector_alloc (3);
		gsl_matrix *evec = gsl_matrix_alloc (3, 3);
		
		gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (3);
 
		gsl_eigen_symmv(&m.matrix, eval, evec, w); 
 
		gsl_eigen_symmv_free(w);
		gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);
		
		
		for (uint i = 0; i < 3; i++) {

			gsl_vector_view evec_i = gsl_matrix_column (evec, i);
			
			if (i==0) {
				C.x = gsl_vector_get(&evec_i.vector, 0);
				C.y = gsl_vector_get(&evec_i.vector, 1);
				C.z = gsl_vector_get(&evec_i.vector, 2);
			}
			
			if (i==1) {
				B.x = gsl_vector_get(&evec_i.vector, 0);
				B.y = gsl_vector_get(&evec_i.vector, 1);
				B.z = gsl_vector_get(&evec_i.vector, 2);
			}
			
			if (i==2) {
				A.x = gsl_vector_get(&evec_i.vector, 0);
				A.y = gsl_vector_get(&evec_i.vector, 1);
				A.z = gsl_vector_get(&evec_i.vector, 2);
			}
			
		}
		
		gsl_vector_free (eval);
		gsl_matrix_free (evec);
		
	}
	
};

#endif